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ABSTRACT 

In  the  computation  of  helicopter  rotor  flow  flelds,  wake  effects  can  be  very 
important  since  each  blade  passes  close  to  the  wake  produced  by  the  preceding 
ones,  causing  a  large  local  effect.  Also,  the  vortical  flow  from  the  wake  of  a  number 
of  blade  passages  causes  a  large  global  effect.  Extensive  work  has  been  done  on 
helicopter  rotor  flows  using  integral  methods  to  follow  the  vortex  wake.  These 
generally  cannot  treat  compressibility  effects.  Also,  they  have  diflBculty  attaining 
stable  solutions,  particulaurly  in  hover,  where  a  large  number  of  interacting  vortex 
sheets  must  be  treated. 

A  new  method  has  been  developed  which  like  integral  methods  does  not  con¬ 
strain  or  spread  the  wake.  Also,  like  finite  difference  methods,  it  can  treat  com¬ 
pressibility  effects.  This  method  has  been  developed  into  a  computer  program  for 
the  computation  of  rotor  flow  fields  in  hover  with  free  wakes.  The  method  utilizes 
a  finite  volume  potential  flow  technique.  The  beisic  appro£w:h  involves  modifying 
the  potential  flow  wake  treatment  so  that,  within  the  niimerical  approximation, 
the  momentiun  is  conserved  there  as  in  the  rest  of  the  flow  field.  The  internal 
structure  of  the  vortex  is  not  solved  for,  but  is  modeled  and  spread  over  several 
grid  points.  The  wake  position  and  vorticity  strength  are  computed  so  that  the 
momentum  over  the  wake  is  balanced  in  an  integral  sense.  Results  computed  by 
this  approach  for  the  circulation  and  wake  geometry  are  compared  with  experi¬ 
mentally  measured  data.  Cases  treated  include  subsonic  and  transonic  flows,  high 
and  low  aspect  ratios,  and  two  and  four-biaded  rotor  configurations. 
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CHAPTER  1 


INTRODUCTION 

The  task  of  predicting  the  flow  field  ajid  airloads  of  a  helicopter  rotor  con¬ 
tinues  to  be  of  primary  importance  for  providing  and  evaluating  improved  rotor 
designs.  Unfortunately,  rotary-wing  flow  fields  are  among  the  most  complex  in 
aerodynamics.  Recently,  application  of  higher  level  computational  aerodynamics 
techniques  to  this  problem  has  become  possible  with  the  advent  of  high  speed  large 
memory  computers.  Computational  studies  have  been  conducted  to  investigate  the 
influence  of  the  rotor  wake  and  develop  methodology  for  predicting  the  rotor  flow 
field  and  blade  airloads  [l].  Most  current  methods  for  computing  rotor  flow  fields, 
including  wakes,  are  based  on  integral  or  finite  difference  techniques  [2j.  Integral 
methods,  in  the  form  of  lifting  line  or  lifting  surface  methods,  are'  able  to  treat  the 
rotor  blades  as  well  as  the  wake  in  a  single  unified  way  3j.  They,  however,  have 
several  disadvantages:  they  are  restricted  to  linear,  low  speed  incompressible  cases 
and  they  have  difficulty  attaining  stable  solutions  in  hover  where  a  large  number 
of  interacting  vortex  sheets  must  be  treated.  Finite  difference  methods  (also  finite 
volume  and  finite  element),  on  the  other  hand,  usually  can  treat  compressibilty 
effects  and  can  accurately  solve  for  flow  in  the  immediate  region  of  the  blade. 
Most  of  these  methods,  which  include  Potential  Flow  and  Euler  based  methods, 
however,  do  not  treat  the  flow  field  in  a  single  unified  way:  they  isolate  the  region 
close  to  the  blade  and  only  solve  the  difference  equations  in  those  regions.  The 
entire  vortex  system  is  currently  treated  in  these  methods  by  coupling  with  an 
integral  technique,  such  as  a  classical  lifting  line  or  lifting  surface  method.  This 
vortex  calculation  is  used  to  determine  an  induced  velocity  on  the  surface  of  the 
finite  difference  grid  surrounding  each  blade,  which  is  then  used  2is  a  boundary 
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condition  for  the  finite  difference  calculation  in  a  coupled  iteration  scheme. 

1.1  Background 

For  hovering  rotor  wake  geometry,  the  fundamental  wake  structure  was  con¬ 
firmed  in  the  1960’s  and  early  1970’s  during  model  rotor  tests  of  rotors  with 
varying  numbers  of  blades,  twist,  taper,  camber,  tip  shape,  thrust  level,  and  tip 
speed  [4].  Smoke  flow  visualization  was  used  to  photograph  cross  sections  of  the 
wake  as  shown  in  the  sample  flow  visualization  photograph  in  Figure  1.  The  fun¬ 
damental  wake  contains  two  primary  components.  The  first,  and  most  prominent, 
is  the  strong  tip  vortex  which  arises  from  the  rapid  rolling  up  of  the  portion  of 
the  vortex  sheet  shed  from  the  tip  region  of  the  blade.  The  second  feature  is 
the  vortex  sheet  shed  from  the  inboard  portion  of  the  blade  (Figure  2).  Flow 
visualization  photographic  data  were  analyzed  to  determine  wake  coordinates  for 
various  rotor  designs  and  operating  conditions.  Generalized  wake  equations  were 
developed  for  the  tip  vortex  and  inboard  wake  geometry,  as  generalized  functions 
of  number  of  blades,  twist,  and  thrust  level  f4,  5j.  The  characteristics  of  the  wake 
geometry  largely  result  from  the  velocities  induced  by  the  strong  tip  vortex.  The 
exact  nature  of  this,  however,  has  been  difficult  to  distinguish  in  flow  visualization 
studies. 

1.2  Rotor  Analysis  in  Hover 

An  understanding  of  rotor  wake  behavior  is  of  primary  importance  in  the 
prediction  of  aerodymamic  loading.  It  is  clear  that  any  complete  analysis  of  the 
rotor  problem  must  account  for  four  distinct  yet  coupled  phenomena: 


Figure  1:  Sample  Flow  Visualization  Photograph  (From  Reference  (41  ) 
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1)  the  nonlinear  compressible  flow  near  the  blade; 

2)  the  highly  complex  vortex  structure; 

3)  close  blade-vortex  interactions; 

4)  viscous  flow  phenomena  such  as  stall. 

Analyses  for  each  of  these  aspects  of  the  problem  already  exist,  though  not  in  a 
rigorously  coupled  form.  A  major  question  in  advanced  aerodynamic  methods  is 
how  can  the  wake  be  included  in  rotor  flow  field  analysis.  A  number  of  approaches 
for  modeling  the  wake  have  been  developed.  In  the  wake  modeling  approach,  the 
analysis  can  be  resolved  into  three  components:  the  aerodynamic  theory  used  for 
the  inflow  solution,  the  method  of  calculating  rotor  performance  from  the  inflow 
solution,  and  the  wake  model. 

By  far  the  most  common  flow  analysis  methods  currently  used  in  the  rotor 
industry  are  various  integral  formulations  of  the  potential  equation.  These  include 
panel  methods,  lifting  line  and  lifting  surface  approaches.  This  class  of  schemes 
provides  the  easiest  and  most  efficient  methods  to  predict  the  inflow  induced  by  the 
rotor’s  complex  wake  system  when  the  geometry  of  the  system  can  be  specified  a 
priori.  The  analysis  based  on  a  lifting  line  representation  for  the  blade  is  reported 
to  predict  hover  performance  for  a  wide  range  of  rotors  !6l.  But  this  representation 
is  not  suitable  for  complex  planforms.  Lifting  surface  methods  have  been  applied 
to  hovering  rotor  analysis  by  a  number  of  investigators  [7,  8].  Typically,  modeling 
of  the  influence  of  the  complete  rotor  and  wake  has  been  accomplished  with  these 
methods.  In  this  approach  the  blade  geometry  is  represented  by  a  thin  surface 
composed  of  a  number  of  conveniently  chosen  trapezoidal  panels,  forming  a  vortex 
circuit  in  the  form  of  a  rectangle.  The  geometry  of  the  wake  is  either  computed 
or  prescribed.  The  wake  is  also  represented  by  a  multiple  vortex  circuit  system. 
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Most  of  the  effects  of  compressibility  camnot  be  obtained  by  this  method. 

Panel  methods  use  a  surface  singularity  distribution  to  solve  the  linear  poten¬ 
tial  equation  for  arbitrary  geometry  [9]  .  The  use  of  a  velocity  potential  reduces 
the  problem  to  the  solution  of  an  integral  equation  for  a  scalar  function.  A  free 
wake  analysis,  which  is  the  most  sophisticated  approach  to  wake  modeling,  con¬ 
siders  the  wake  as  a  free  vortex  flow.  A  number  of  the  panel  method  solutions  for 
the  hovering  rotor  include  a  free  wake  geometry  calculation  [10,11,12].  In  these, 
the  wake  is  divided  into  near,  intermediate,  and  far  wake  regions.  Vortex  filaments 
are  distributed  in  the  near  and  intermediate  wake  regions  and  an  analytical  model 
used  to  represent  the  far- wake. 

« 

Miller  [11]  developed  a  hover  free  wake  geometry  calculation,  using  two- 
dimensional  and  axisymmetric  models  for  efficiency.  The  wake  was  assumed  to 
have  rolled  up  into  three  vortices  (tip,  inboard,  and  root)  which  were  replaced 
by  a  far  wake  model  after  four  revolutions.  Murman  and  Stremel  12]  calculated 
two-dimensional  unsteady  wake  development  by  a  cloud-in-cell  method.  Discrete 
vortex  elements  in  the  wake  were  tracked  in  a  Lagrangian  frame.  For  each  time 
step,  the  vorticity  was  distributed  to  a  fixed  mesh,  on  which  the  velocities  were 
calculated  by  a  finite  difference  solution  of  Laplace’s  equation.  Stremel  [13]  applied 
this  two-dimensional  model  to  hovering  rotor  wake  roll-up.  Roberts  [I4j  developed 
a  cloud-in-cell  method  to  calculate  hovering  rotor  wake  geometry. 

In  recent  years  a  number  of  finite-difference  codes  for  the  prediction  of  tran¬ 
sonic  rotor  flows  has  been  developed.  Egolf  and  Sparks  [15]  coupled  a  full  potential 
solution  for  the  near  field  of  a  hovering  rotor  blade  with  an  incompressible  discrete 
vortex  solution  for  the  far  field.  The  wake  geometry  waa  prescribed,  and  the  far 


7 


field  solution  determined  the  velocity  potential  on  the  boundaries  of  a  finite  differ¬ 
ence  computational  domain.  Sauikar  et  al  [16]  solved  three  dimensional  unsteady 
Euler  equations  for  a  hovering  rotor  in  a  small  box  around  each  blade.  Caradonna 
and  Strawn  [17]  solved  the  full  potential  equation  in  conservative  form  also  in  a 
box  around  each  blade  and  Incorporated  a  wake  modeling  scheme.  The  basic  ap¬ 
proach  to  this  wake  modeling  was  to  separate  the  rotor  wake  system  into  two  parts. 
The  most  important  part,  the  tip  vortices,  were  explicitly  modeled  in  the  finite 
difference  solution  procedure  when  they  were  inside  the  computational  box.  The 
second  part  of  the  wake  system,  consisting  of  vorticity  elements  that  were  outside 
the  computational  grid,  was  modeled  with  an  inflow  boundary  condition  at  the 
boundary  of  the  computational  grid.  A  similar  method  was  developed  by  Tung 
and  Chang  [18|  with  a  non-conservative  formulation  of  the  potential  equation. 

The  approach  used  in  these  finite-difference  methods  to  include  the  total  wake 
influence  has  been  to  include  the  effects  of  near  vortices  in  the  form  of  an  inflow 
boundary  condition  applied  on  the  surface  [l7]  with  the  inflow  being  determined  by 
a  Biot-Savart  computation.  The  outer  vorticity  contributions  are  included  either 
on  the  grid  outer  boundary  or  specified  as  a  spanwise  varying  correction  to  the 
blade  twist  distribution. 

The  investigation  presented  herein  describes  a  unified  method  that  is  fully 
compressible  and  which  computes  the  wake  without  requiring  external  specification 
of  the  wake,  or  separate  computations  for  the  wake  and  blade  region.  The  problems 
associated  with  other  methods  are  avoided  by  using  an  embedding  procedure  by 
which  vorticity  layers  can  be  put  anywhere  in  the  grid.  The  method  is  based 
on  the  fact  that  any  flow  field  can  be  decomposed  into  potential  and  vortical 
parts.  A  potential  is  defined  on  a  set  of  grid  points,  as  in  standard  potential 
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methods.  Also,  as  in  these  methods,  mass  balance  relations  are  enforced  at  each 
point  on  the  grid.  The  potential,  however,  does  not  have  any  discontinuities, 
and  therefore  does  not  represent  any  vorticity.  The  vorticity  is  represented  by  a 
separate  velocity  field  which  is  added  to  the  gradient  of  the  potential.  The  location 
and  strength  of  this  velocity  field  is  determined  by  momentum  considerations. 
This  added  field  is  concentrated  in  sheets,  as  is  the  vorticity  that  results  from 
it.  No  external  specification  of  the  strength  of  these  sheets  or  their  location  is 
required:  momentum  conservation  relations  together  with  mass  balance  relations 
are  used  to  determine  these.  This  is  believed  to  be  the  first  complete  treatment  of 
this  problem.  The  method  utilizes  a  modification  of  a  compressible  finite  volume 
potential  flow  technique  frequently  used  in  fixed-wing  analysis,  which  has  been 
developed  into  a  cpmputer  program  for  the  computation  of  compressible  rotor 
flow  fields  in  hover  with  free  wakes.  Free  wake  geometry  determined  with  this  new 
wake  treatment  are  in  substantial  agreement  with  experiment. 
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CHAPTER  2 


GENERAL  ANALYSIS 


2.1  Elementary  Considerations  of  Hovering  Rotor 


In  the  momentum  theory  analysis  the  rotor  is  modeled  as  an  actuator  disk, 
which  is  a  circular  surface  of  zero  thickness  that  can  support  a  pressure  difference. 
The  flow  problem  is  solved  by  using  the  basic  conservation  laws  of  fluid  motion. 
The  actuator  disk  model  is  only  an  approximation  to  the  actual  rotor.  Momentum 
theory  gives  the  induced  power  per  unit  thrust  for  a  hovering  rotor: 

/  T 

V^nduced  =  ^/^  =  y  ^ 

This  relation  determines  the  basic  characteristic  of  helicopter  rotor.  Simple  mo¬ 
mentum  theory  provides  some  understanding  of  the  basic  relationships  between 
such  important  design  parameters  as  disc  loading  and  power  required  per  unit  of 
thrust.  However,  this  theory  has  serious  limitations  in  providing  guidance  for  rotor 
design  and  does  not  provide  a  physical  concept  that  could  explain  the  nonunifor¬ 
mities  of  downwash  velocities. 


Vortex  theory  uses  the  Biot-Savart  law  for  the  velocity  induced  by  the  wake 
vorticity.  In  this  approach  the  rotor  is  modeled  by  segments  of  a  vortex  filament. 
The  vortex  elements  are  allowed  to  be  transported  by  the  resultant  velocity  of 
the  free  stream  and  vortex  induced  velocities  calculated  using  the  Biot-Savart  law. 
The  computational  procedure  required  in  the  free  wake  analysis  is  huge  and  can 
be  accomplished  only  by  the  use  of  high  capacity  computers.  In  addition,  various 
safeguards  must  be  used  to  avoid  mathematical  singularities  that  have  no  physical 
counterpart. 
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In  recent  years  a  substantial  number  of  finite-difference  codes  for  the  predic¬ 
tion  of  three-dimensional  transonic  flows  have  been  developed.  These  codes  have 
been  successfully  used  in  fixed-wing  flow  analysis.  The  development  and  use  of 
corresponding  rotary-wing  codes  has  proceeded  very  slowly.  A  primary  reason  for 
this  delay  is  the  much  greater  complexity  of  the  rotor  environment. 

Classical  lifting-line  theory  treats  the  case  of  a  high  aspect-ratio,  planar,  fixed 
wing  in  steady  flow.  In  the  linearized  model  both  the  wing  and  wake  are  repre¬ 
sented  by  this  planar  sheet  of  vorticity.  For  the  fixed  wing  the  distortion  of  the 
wake  and  the  rollup  of  tip  vortices  can  generally  be  neglected,  because  the  wake  is 
convected  downstream  away  from  the  wing.  The  basic  assumption  of  lifting-line 
theory,  that  the  wing  has  high  aspect  ratio,  can  sometimes  be  satisfied  with  a 
rotor  blade.  However,  the  requirement  that  there  be  no  rapid  change  in  spanwise 
circulation  restricts  the  validity  of  lifting-line  theory  rotor  flow  even  though  the 
geometric  aspect  ratio  is  large.  There  are  two  important  regions  of  general  rotor 
flow  fields  where  this  requirement  is  not  satisfied:  at  the  blade  tip  and  near  a 
blade-vortex  interaction.  On  rotor  blades  the  loading  is  typically  concentrated 
near  the  tip  because  of  the  rotation  of  the  blade,  and  the  gradient  of  the  lift  is 
particularly  high  there.  Hence,  any  small  distortion  of  the  loading  due  to  the  three 
dimensional  flow  effect  is  very  important.  Conservation  of  vorticity  requires  that 
the  bound  circulation  be  trailed  into  the  wake  from  the  blade  tip  and  root.  The 
vorticity  is  also  left  in  the  rotor  wake  as  a  consequence  of  radial  and  azimuthal 
change  in  the  bound  circulation  (Figure  3).  In  classical  theory,  the  trailing  vor¬ 
ticity,  Tft,  which  is  due  to  the  radial  variation  in  bound  circulation,  is  parallel  to 
the  free  stream.  The  shed  vorticity,  'jg,  due  to  the  azimuthal  variation  in  bound 
circulation,  is  oriented  radially  in  the  wake.  The  strength  of  the  rotor  trailing  and 


Figure  3;  Schematic  of  Trailed  and  Shed  Vorticity  in  Rotor  Wake  (From 
Reference  [191  ) 


12 


shed  vorticity  we  given  by, 


Is  =  - 


J_^ 


Hover  is  the  operating  state  in  which  the  lifting  rotor  has  no  forward  velocity 
relative  to  the  air.  This  implies  axial  symmetry  of  the  rotor,  and  hence,  that 
loads  on  the  rotor  blades  are  independent  of  the  azimuth  position.  Hence,  this 
particular  shed  vorticity  is  not  of  importance  in  hover  airload  computations.  While 
the  circulation  drops  to  zero  over  a  finite  distance,  the  rate  of  decrease  is  still  very 
high.  The  result  is  a  large  trailing  vorticity  strength  at  the  outer  edge  of  the  wake, 
causing  the  vortex  sheet  to  roll  up  into  a  concentrated  tip  vortex.  The  tip  vortex 
has  a  small  core  radius  that  depends  on  the  blade  geometry  and  loading.  On  the 
inboard  portion  of  the  blade,  the  bound  circulation  drops  off  gradually  to  zero  at 
the  root.  Hence  there  is  an  inboard  sheet  of  trailed  vorticity  in  the  wake  with 
opposite  sign  to  the  tip  vortex.  Since  the  gradient  of  the  bound  circulation  is  low, 
the  root  vortex  is  generally  much  weaker  amd  more  diffused  than  the  tip  vortex. 
Hence,  the  strong  concentrated  tip  vortices  are  by  fax  the  dominant  feature  of 
the  rotor  wake.  In  addition,  because  of  its  rotation,  a  rotor  blade  encounters  the 
tip  vortex  from  the  preceding  blade.  When  a  vortex  passes  close  to  the  blade  it 
induces  a  large  velocity  and  hence  a  large  change  in  loading  on  the  rotor.  Methods 
have  been  developed  to  calculate  the  nonuniform  intlow  due  to  the  wake,  the  self 
induced  distortion  of  the  wadce,  and  the  vortex  induced  loads  flQ].  From  the 
dominant  role  of  the  tip  vortices  in  determining  the  inflow  and  loading,  it  follows 
that  the  determination  of  its  position  is  the  most  important  part  of  the  rotor  wake 
geometry  calculation. 
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2.2  Hovering  Rotor  Wake  Models 

There  are  a  number  of  methods  used  to  specify  wake  geometry  in  rotary 
wing  analysis.  The  rigid  wake  model  implies  an  undisturbed  helical  geometry,  in 
which  all  the  wake  elements  are  convected  with  the  same  mean  axial  velocity.  In 
this  model,  the  wake  consists  of  helical  vortex  sheets  that  move  as  rigid  surfaces 
with  uniform  velocity  and  no  distortion.  The  transport  velocity  of  the  wake  is 
determined  by  rotor  disk  loading,  and  the  pitch  of  the  helical  surface  is  determined 
by  the  axial  and  rotational  velocity  of  the  blades.  The  wake  movement  imparts 
a  velocity  to  the  fluid  at  the  wake  surface.  In  the  limit  of  an  infinite  number  of 
blades,  the  sheets  are  very  close  and  as  a  consequence  all  the  fluid  is  carried  with 
the  wake  and  there  are  no  losses  due  to  flow  around  the  edges.  With  a  finite 
number  of  blades  there  will  be  radial  flow  as  well,  which  decreases  the  lift  at  the 
blade  tip.  An  elementary  extension  of  this  model  is  the  semi-rigid  wake  model,  in. 
which  each  element  is  convected  downward  with  the  induced  velocity  at  the  point 
on  the  rotor  where  it  was  created.  The  simplest  wake  model  including  the  effects  of 
a  finite  number  of  blades  consists  of  helical  vortex  sheets  trailed  from  each  blade. 
The  major  effect  of  the  finite  number  of  blades  is  a  reduction  of  the  loading  at  the 
blade  tip.  This  effect  is  generally  represented  by  a  tip  loss  correction  for  a  finite 
number  of  blades. 

A  wake  model  consisting  of  undistorted  vortex  sheets  is  not  adequate  for  he¬ 
licopter  rotors,  where  the  interference  between  the  rotor  blade  and  wake  vorticity. 
and  the  self  induced  distortion  in  the  wake,  are  important.  The  trailed  vortic¬ 
ity  quickly  rolls  up  into  concentrated  tip  vortices  that  remain  close  to  the  rotor 
and  strongly  influence  the  loading  near  the  tips  of  both  generating  and  following 
blades.  The  most  accurate  specification  of  the  wake  is  the  free  wake  model.  This 


includes  the  distortion  from  the  basic  helix,  as  each  wake  element  is  convected 
with  the  local  flow,  including  the  velocity  induced  by  the  wake  itself.  When  mea¬ 
sured  wake  geometry  information  is  used,  instead  of  a  calculation  this  is  called  a 
prescribed  wake  model.  Clark  and  Leiper  [20]  developed  a  method  for  calculating 
the  distorted  wake  geometry  of  hovering  rotors.  Their  wake  model  consisted  of  a 
number  of  constant  strength  trailed-vorlex  lines.  The  far  wake  was  approximated 
by  segments  of  ring  vortices.  Two  revolutions  of  the  free  wake  were  used,  followed 
by  30  revolutions  of  far  wake.  Landgrebe  [4j  conducted  an  experimental  investiga¬ 
tion  of  the  performance  and  wake  geometry  of  a  model  hovering  rotor.  The  wake 
geometry  was  measured  by  flow  visualization,  and  the  data  was  used  to  develop 
expressions  for  the  axial  convection  and  radial  contraction  of  the  tip  vortices.  The 
tip  vortex  elements  were  found  to  have  a  roughly  constant  descent  rate  before  and 
after  passing  beneath  the  following  blade.  Prior  to  the  encounter  with  the  blade 
the  descent  rate  is  approximately  proportional  to  the  blade  loading,  C^/cr.  After 
the  encounter  the  axial  convection  rate  is  higher  and  is  approximately  propor¬ 
tional  to  the  mean  inflow,  This  generalized  wake  geometry  information 

was  used  in  rotor  performance  calculations.  Scully  [21]  developed  a  free  wake  anal¬ 
ysis  method  for  a  helicopter  rotor.  The  procedure  for  calculating  wake  geometry 
consisted  basically  of  integrating  the  induced  velocity  at  each  wake  element.  Once 
the  wake  geometry  is  specified,  either  a  lifting  surface  or  a  lifting-line  theory  can 
be  used  to  obtain  blade  loading. 

2.3  Potential  Flow  Methods 

Lifting  surface  theory  retains  the  mutual  interaction  of  ail  elements  of  the 
rotor  and  wake  by  representing  the  rotary-wing  by  vortex  surfaces  and  by  satisfying 
appropriate  boundary  conditions  over  the  entire  surface.  A  free  wake  geometry 
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calculation  is  required,  and  also  a  calculation  of  the  tip  vortex  roll-up  and  other  fine 
structures  of  the  wake.  Hence,  the  computation  time  required  increases  rapidly 
with  the  number  of  elements  used  to  represent  the  blade/vortex  system. 

From  the  above  discussion  of  these  vortex  theories,  it  can  be  seen  that,  though 
they  offer  a  precise  description  of  the  rotor  aerodynamics  in  practical  application, 
they  require  a  large  amount  of  computational  effort.  The  application  of  a  velocity 
potential  which  will  be  described  below,  makes  it  possible  to  determine  steady  and 
unsteady  flow  fields  induced  by  the  rotor  in  both  incompressible  and  compressible 
fluids  with  a  precision  similar  to  that  offered  by  the  above  vortex  theories  but  with 
less  computational  effort. 

2.3.1  Velocity  Potential  in  Incompressible  Flow 


The  continuity  equation  for  an  incompressible  fluid  requires  that 

du  dv  dw 

- 1 - i - =  0 

dx  dy  dz 

further  restricting  this  perfect  fluid  to  be  irrotational  leads  to  the  existence  of  a 
velocity  potential,  <t>.  The  components  of  the  velocity  vector  associated  with  4>  can 
be  obtained  as 

v(F)  =  •V0(F) 


The  condition  of  continuity  can  then  be  written  as 


V  • 


_  d^4>  d^<i> 

dx^  dy^  ^  dz"^ 


resulting  in  Laplace’s  equation  for  <i>. 
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This  equation  relates  the  potential  at  an  arbitrary  point  in  the  fluid  to  the 
superposition  of  potentials  due  to  doublet  and  source  singuluities  distributed  on 
the  body  surface.  Kocurek  [21)  applied  this  singularity  type  of  discretization  to 
solve  the  problem  of  a  hovering  rotor. 

2.3.2  Potential  Methods  in  Compressible  Fluids 

At  present  a  number  of  schemes  axe  being  developed  to  compute  compressible 
flows.  A  variety  of  approaches  use  potential  methods.  The  transonic  flow  past  a 
non-lifting,  hovering  rectangular  rotor  was  first  solved  by  Caradonna  and  Isom 
[22).  They  developed  a  finite  difference  calculation  procedure  for  the  Transonic 
Small  Disturbance  equation.  Later,  it  was  extended  to  solve  the  Full  Potential 
equation.  For  lifting  rotors,  the  conventional  potential  flow  method  was  not  able 
to  treat  the  rotor  blade  m  well  as  the  wake  in  a  single  unified  way.  In  this  method, 
the  vortex  sheets  are  treated  as  discontinuities  in  the  potential,  which  constrains 
the  sheets  to  lie  on  segments  of  surfaces  of  the  computational  grid.  Accordingly,  in 
conventional  potential  methods,  only  short  segments  of  the  sheet,  which  are  fairly 
flat,  and  where  the  grid  can  be  distorted  to  follow  the  sheet,  caji  be  accurately 
treated. 

More  recently,  a  coupled  iteration  scheme  has  been  developed  where  a  finite 
difference  solution  procedure  is  obtained  in  a  small  region  near  each  blade.  The 
computation  In  this  grid  is  coupled  to  an  exterior  vortex  filament  wake  model, 
which  describes  almost  all  of  the  wake  influence.  Unfortunately  there  are  major 
drawbacks  to  this  approach.  First,  unlike  the  computation,  there  is  no  clear  sep¬ 
aration  in  the  physical  flow  between  the  wake  and  the  flow  in  the  blade  region.  If 
accurate  solutions  are  to  be  obtained,  the  flow  region  that  is  computed  using  the 
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compressible  method  on  the  grid  must  extend  a  large  number  of  chords  from  the 
blade.  The  wake  then  typically  penetrates  this  region  and  must  be  treated  within 
the  grid.  Since  the  size  of  this  region  as  well  as  the  wake  geometry  near  the  blade 
will  depend  sensitively  on  the  flow  conditions,  deciding  which  portion  of  the  wake 
spirals  must  be  included  in  the  grid  can  be  complicated  for  the  general  cases.  Fur¬ 
ther,  simply  embedding  the  vortex  sheet  as  a  potential  discontinuity  or  embedded 
line  vortices  is  only  feasible  for  simple  cases,  such  as  where  the  sheet  rolls  up  into 
a  single  well-defined  vortex  near  the  tip  and  one  near  the  root.  Further  drawbacks 
concern  the  feasibility  of  treating  the  general  wake,  even  outside  of  the  grid  re¬ 
gion,  as  a  collection  of  vortex  filaments.  While  this  is  possible  in  simple  cases,  in 
general  it  is  very  difficult  since  the  vortical  regions  can  leave  the  blade  with  even 
the  signs  of  vorticity  varying  along  the  blade.  Also  if  Jthere  are  a  large  number  of 
filaments,  these  must  be-integrated  over  each  time  step  to  determine  the  velocity 
induced  by  the  wake  over  the  outer  boundaries  of  the  grid  regions.  Accurately 
interfacing  this  calculation  with  a  set  of  grids  near  each  blade  would  appear  to 
be  complicated.  Furthermore,  a  filament  approach  complicates  the  computation 
of  near  blade-vortex  interactions  because  special  smoothing  must  be  applied  near 
each  filament  to  eliminate  locally  infinite  velocities. 

2.3.3  Vortex  Embedding  Method 

The  basic  approach  described  herein  involves  modifying  the  potential  flow 
wake  treatment  so  that,  within  the  numerical  approximation,  momentum  is  con¬ 
served  there,  as  in  the  rest  of  the  field.  The  internal  structure  of  the  vortex  sheet 
is  not  solved  for,  but  is  modeled  and  spread  over  several  grid  points.  The  wake 
position  and  vorticity  strength  are  computed  so  that  momentum  over  the  wake  is 
balanced  in  an  integral  sense.  The  potential  is  defined  on  a  set  of  grid  points  as  in 
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standard  methods.  Also,  as  in  these  methods,  mass  balance  relations  are  enforced 
at  each  point  on  the  grid.  However,  the  potential  does  not  have  any  discontinu¬ 
ities,  and  therefore  does  not  represent  any  vorticity.  The  vorticity  is  represented 
by  a  separate  velocity  field  which  is  added  to  the  gradient  of  the  potential.  This 
added  field  is  concentrated  in  sheets,  as  is  the  vorticity  that  results  from  it.  No 
external  specification  of  the  strength  of  these  sheets  or  their  locations  is  required: 
momentum  conservation  relations  are  used  to  determine  these,  together  with  the 
mass  balance  relation.  Thus,  it  is  a  true  free-wake  method. 


Initial  use  of  the  method  for  a  vortex  line  conyecting  past  a  wing  in  compress¬ 
ible  flow  was  published  in  1983  [22].  A  preliminary  use  of  the  method  for  roll-up 
of  vortex  sheets  was  also  published  in  1983  [23] . 


2.3.4  Description  of  the  Vortex  Embedding  Procedure 

The  basis  of  the  Vortex  Embedding  procedure  lies  in  the  fact  that  any  velocity 
q  can  be  expressed  as  the  sum  of  potential  and  vortical  components.  First  the 
velocity  is  decomposed  into  a  free  stream,  potential,  and  vortical  part: 


9  = 


=  w  X  f  +  V(^  + 


(1) 


where  w  is  the  rotor  angular  velocity,  directed  along  the  axis,  f  is  normal  to  the 
axis  and 

u;  =  V  X  9 


The  vortical  part,  9  ",  is  concentrated  near  the  sheet.  This  decomposition  is  made 
especially  powerful  by  the  fact  that  for  a  given  uJ  distribution  7  ”  is  not  unique. 
This  can  produce  some  very  useful  simplifications  in  various  flow  models.  For 
instance,  in  a  typical  potential  computation  for  the  lifting  flow  about  a  wing,  the 
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entire  velocity  is  usually  represented  by  the  potential,  (that  is  ^  "  =0)-  However, 
the  production  of  lift  implies  the  shedding  of  vorticity.  With  a  potential  alone  this 
can  only  be  represented  by  a  sheet  along  which  <t>  is  discontinuous  and  where  the 
magnitude  of  this  discontinuity  is  F,  the  shed  circulation.  This  sheet  typically  must 
lie  along  a  grid  coordinate.  However,  the  use  of  velocity  decomposition  allows  a 
considerable  simplification  and  versatility  in  the  sheet  description.  The  basic  idea 
is  to  spread  the  sheet  into  a  thin  vorticity  region  described  by  an  appropriate 
velocity  distribution,  g 

A  fixed  grid  (in  the  rotating  blade-fixed  frame)  is  used  to  solve  the  compress¬ 
ible-potential  equation  for  the  potential,  <f>: 


dcifiu)  +  dy(pv)  -h  d,{pw)  =  0  (2) 

where  p  is  the  density  and  has  the  isentropic  form  away  from  the  sheet  based  on 
the  total  velocity,  q,  with  components  u,  v  and  w. 


P  = 


(3) 


The  vortical  component,  q'’,  is  spread  over  several  grid  points  around  the  vortex 
sheet  so  that  the  vorticity  is  concentrated  there,  q  "  satisfies 


Vxg" =w 


During  the  iteration  towards  convergence  (the  solution  is  steady  in  the  blade 
coordinate  system  for  hover)  a  four  step  procedure  is  used  repeatedly: 


1.  The  vortex  position  is  integrated  as  a  set  of  marker  streamlines  to  follow 
the  flow  using  interpolated  value  of  q  from  fixed  grid; 


20 


2.  ^  "  is  computed  at  grid  points  near  the  sheet; 

3.  A  potential,  <l>  is  computed  at  each  grid  point  to  satisfy  Equation  (2) 
with  q  "  fixed. 

4.  A  new  velocity  q\s  computed  at  each  grid  point  based  on  Equation  (1). 

At  convergence  Equation  (2)  is  satisfied  and  the  vortex  sheet  follows  the  flow. 
Although  the  vortex  sheet  is  effectively  spread  over  several  grid  points,  q"  is 
a  continuous  vortical  velocity  field  which  suffers  no  artificial  numerical  diffusion. 
The  internal  structure  of  the  sheet  can  be  modeled  by  simply  choosing  appropriate 
values  or  solved  for  using  a  viscous  flow  solution.  Also,  experimental  or  theoretical 
information  can,  if  desired,  be  inserted  into  the  method  to  prescribe  this  internal 
structure,  while  no  experimental  information  is  needed  for  the  wake  position. 
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CHAPTER  3 
SOLUTION  PROCEDURE 


3.1  Basic  Formulation 

A  potential  flow  solver  using  a  rapidly  converging  approximate  factorization 
scheme,  together  with  a  finite  volume  discretization  scheme  is  employed.  The 
scheme,  originally  developed  for  a  fixed  wing,  is  described  in  Reference  [24,  25j. 

First,  a  mapping  is  defined  from  physical  (x,  y,  z)  to  computational  (X,  Y,  Z) 
coordinates.  A  set  of  constant  radial  computational  planes  are  defined.  In  each 
plane  the  airfoil  is  mapped  to  a  segment  on  the  real  axis  using  the  transformation 
described  in  Reference  [26j.  If  the  transformations  are  represented  by  a  Jacobian 
matrix  H,  Equation  (3)  may  be  written  in  terms  of  transformed  variables: 

=  0  (4) 

where  U,  V,W  are  contravariaint  velocities,  and  h  is  the  transformation  Jacobian, 


h  =  det  (H) 


The  physical  velocity  components,  in  terms  of  the  potential  and  vortical  velocities, 
are: 

(5) 


q=\  ^  j  ~  |^(^V.j-rnxf+7'' 

where  only  derivatives  of  the  potential  with  respect  to  the  computational  coordi¬ 
nates  are  required.  The  contravariant  velocities  aire  then  calculated  in  terms  of 
physical  velocities. 


H  is  given  by 


/  A 

^  =  I  tfx  Vy  Vz 
V  ^x  ^z  } 

Equation  (4)  can  now  be  represented  as  a  flux  balance  relation.  For  this  purpose,  a 
secondary  set  of  cells  is  introduced  interlocking  with  the  primary  cells  (Figure  4). 
In  the  computational  domain  the  faces  of  the  secondary  cells  span  the  midpoints 
of  the  primary  cells.  The  purpose  of  the  secondary  cells  is  to  serve  as  control 
volumes  for  the  flux  balance.  The  formula  for  the  local  flux  balance  can  now  be 
written  using  a  second  application  of  the  scheme  on  the  secondary  cells,  where 
Equation  (4)  is  approximated  by 

^YZ<^x(<»^^)  +/^ZX<5v(MI^)  +  =0  (6) 

where  the  quantities  phU,  phV ,  and  pkW  are  the  fluxes  across  the  faces  of  the 
secondary  cell.  This  forrhuta  is  equivalent  to  calculating  the  flux  across  the  part 
of  a  face  of  a  secondary  cell  lying  in  a  particular  primary  cell  using  the  values  for 
p,h,U,V,W  calculated  at  the  center  of  that  primary  cell. 

3.2  Grid  Generation 

The  grid  is  a  particular  type  of  “H”  mesh.  A  two-dimensional  grid  is  generated 
for  each  radial  plane.  The  computational  plane  is  first  covered  by  an  equally-spaced 
grid  in  the  region 

|X|  <  1-5 

|y|<i-5 


where  5  >  0  is  a  small  parameter  that  controls  the  far-field  boundary  position.  A 


stretching  is  performed  that  maps  X  and  Y  from  3:1  to  ±  oo. 


where. 


l+c/2  <  X  <  I  -  6) 
X'  =  X,  [-c/2  <  X  <  c/2) 

Y'  =  1/2K(1  -  Y^) 

X*  =  X  r  c/2 


and  the  airfoil  is  mapped  to  a  segment  on  the  real  axis: 


r=0  [-c/2<X<c2, 


A  square  root  transformation 


5  =  z'-’  =  (X'  ^ 


is  used  to  map  the  Z  plane  to  the  upper  half  plane  (Z). 


The  same  transformation  is  then  used  to  map  the  airfoil  to  a  curve  about 
the  real  axis  in  the  z  plane.  Next,  a  shearing  transformation,  f[z),  is  defined 
in  the  lTn{z)  direction  such  that  this  is  mapped  to  a  segment  of  the  real  axis. 
The  physical  plane  coordinates  are  then  computed  by  subtracting  f(z)  from  the 
z  plane  coordinates  of  the  computational  points.  To  improve  the  resolution  near 
the  leading  edge,  a  nonorthogonal  transformation  is  used  in  the  Z  plane  in  place 
of  the  basic  square  root  one: 


R'  a  yjc  ,2 
=  - =3?  - —  e  ' 

a  ~  \  R  \  ^ 


where  a  is  a  grid  compression  length  scale. 
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This  grid  is  then  rotated  to  the  desired  pitch  angle  as  described  by  the  twist 
distribution  of  the  rotor  blade.  Finally,  a  simple  blending  technique  [27j  is  used 
to  enforce  the  outer  and  periodic  boundary  conditions  on  the  grid  (Figure  5). 

3.3  Discretization 

The  finite  volume  scheme  of  Reference  [28j  is  used  to  discretize  Equation  (6). 
Two  staggered  grids  are  used.  On  one,  p,u,v,w,U,V,W,  H,&nd  h  are  defined;  on 
the  other,  i, y,  z,  X,  K,  Z,  and  Lo{4>)  are  defined.  The  nodes  of  each  grid  are 
at  the  centers  of  the  cell  of  the  other,  and  a  box  scheme  is  used  to  compute  the 
derivatives  in  Equation  (5)  and  Equation  (6).  For  a  variable  /,  we  define 

~  ~  ftj.k+i.  ~  fi,j  +  l,k  ~  fi,3,k\j^ 

Similar  expressions  can  be  written  for  the  Y  and  Z  derivatives. 

The  basic  finite  volume  scheme  leads  to  an  odd-even  decoupling  of  solutions. 
If  a  small  (higher  order)  term  is  added  to  the  right  hand  side  of  the  Equation  (6), 
this  problem  can  be  eliminated.  This  can  be  easily  seen  by  considering  the  case  of 
incompressible  flow  in  two  dimensions.  Setting  h  =  l,p  =  1,  Equation  (6)  reduces 
to 

Mvr  r  =  0 

where  odd  and  even  points  are  decoupled.  It  is  due  to  the  evaluation  of  flux  across 
the  face  labeled  AB  in  Figure  6  using  a  value  of  <t>.^  calculated  at  the  point  A. 
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Figure  5:  Blended  Rotor  Grid 
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If  a  compensation  flux,  €^Y(i)^y,  is  added  across  AB,  the  point  at  which  is 
effectively  evaluated  is  shifted  from  A  to  B.  The  addition  of  similar  compensation 
terms  on  all  faces  produces  the  scheme 

^xx  <f>  f^xx  ^rr  ~  ^^xxrr  ^  =  0 

The  compensation  terms  for  Equation  (6)  can  be  found  in  a  similar  manner  [28]. 
These  are 

B  =  -  V^/a^)  (7) 

C  -W\/a^) 

Where  g,.,  are  the  diagonal  elements  of  {H^ 

Defining  the  following  quantities 

Qxr  =  (>4  -t-  B)n^6yy<p 
Qyz  ={B^C)pL^6y,<i> 

QzX  =  (<^  +  A)Hy6^^(i> 

Qxyz  =  {a  +  b  +  c)6^y^(i> 

Equation  (6)  can  be  written  with  recoupling  terms  as: 

fiy^S^iphU)  +  pg^6y{phV)  +  n^y6^(phW)  + 

^[t^Z^XYQxY  ~  f^X^YZ^YZ  '^f^Y^ZxQzX  ~~  2  ^XYzQxYZl  (®) 

where 

^XYzf  ~  +  ~  fi,]  +  l,k  +  \  ~  fi-^-l,},k-^l,  fi,j,k  +  \ 


fi-i-l,]-rl,k  +  fx,j  +  l,k  +  /t  +  l,;,fc  /t.j.fcj 
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0  <  £  <  1/2 

For  stability  in  supersonic  zones  either  a  first  order  or  a  second  order  artificial 
viscosity  can  be  used.  For  the  first  order  form  a  switching  function  is  defined: 

s  =  max  (0,1  -  -^) 

where  Me  is  the  switching  Mach  number,  taken  to  be  slightly  less  than  1. 

Then,  from  Reference  [25l, 


r2x2. 


uv , 

WU , 

4 

vw 

UV  , 

4 

wu , 

VW  .  . 

(5-  y  0  + 

4 

/X  =  eshp^  (a^ 

In  this,  P,  Qi  and  R  are  constructed  so  that 


P  approximates  —  p,\U' 6^ p 


Q  approximates  -  p.\V\SyP 
R  approximates  -  p,\W\6^p 

In  the  numerical  scheme  Equation  (8)  is  modified  by  the  addition  of  the  terms 


6^P  ■‘r  ByQ  +  6 ^R 

Finally,  an  iterative  scheme  has  to  be  devised  for  solving  the  discretized  equation. 
This  is  accomplished  by  embedding  the  steady  state  equation  in  an  artificial  time 
dependent  equation  [29].  Such  a  process  has  to  be  carefully  controlled  to  ensure 
stability  and  convergence.  It  is  helpful  to  regard  the  iterations  as  successive  levels 
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in  an  artificial  time  coordinate.  The  essential  idea  is  that  the  iterative  procedure, 
regarded  as  a  finite  difference  scheme  for  a  time  dependent  process,  should  be 
consistent  with  a  properly  posed  initial  value  problem  which  is  ‘compatible’  with 
the  steady  state  equation;  that  is  to  say,  that  the  solution  of  the  initial  value 
problem  should  reach  an  equilibrium  point  which  represents  a  solution  of  the 
steady  state  equation.  Then,  if  the  difference  scheme  is  also  stable,  it  should 
converge  to  the  desired  solution. 


Considering  first  the  two  dimensional  case,  the  full  potential  equation  can  be 
written  in  canonical  form: 

~  M^)il>ss  +  4>nn  ^  0  (9) 

where  s  and  n  are  coordinates  in  the  local  stream  and  normal  directions  and 


=■  +2uv<i>xy  + 

0nn  =  \{v^<i>xx  -  2uv4>xy  -r  U^<t>yy) 


In  the  finite  difference  form  for  Equation  (9),  let  updated  values  of  potential 
4>  at  any  level  of  the  calculation  be  denoted  by  the  superscript  Then  the  typical 
central  difference  formula  for  a  second  derivative  in  the  streamwise  direction  can 
be  written  as 


<t>zx  = 


-  ~  +  rAx)(l>:^j  -  (1  -  rAx)d>.y  +  di+i.; 
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Where  the  updated  value  is  used  on  one  side,  the  old  values  on  the  other  side 
(because  the  updated  value  is  not  yet  available),  and  a  linear  combination  of  the 
two  is  used  at  the  center  point.  In  the  time  dependent  system  this  formula  may 
be  interpreted  as  representing: 

At,. 

<Pzx  —  ■ 


With  a  similar  representation,  Equation  (9)  can  be  recast  as  a  time  dependent 
equation  which  contains  mixed  space-time  derivatives; 

(1  -  +  <f>nn  -  2l34>nt  =  0 

The  values  of  a  and  /3  depend  on  the  combination  of  new  and  old  values  actually 
used  in  the  difference  formulas.  While  interpreting  the  difference  scheme  as  the 
representation  of  a  time  dependent  process,  it  can  be  seen  from  stability  analysis 
that  the  presence  of  a  (t>t  term  can  lead  to  instability  at  supersonic  points  [30|. 
This  is  in  contrast  to  the  subsonic  part  of  the  flow.  There,  the  damping  due  to  <i>t 
plays  a  critical  rule  in  the  convergence  of  the  relaxation  scheme. 

The  relaxation  schemes  have  to  be  devised  based  on  these  principles,  extended 
to  three  dimensions.  For  our  equations,  we  solve  a  discrete  approximation  to 

6^  {phU  -hP)+  6y  {phV  +  Q)  +  6^  {phW  +  R)  =  (9) 

In  this  care  should  be  taken  to  have  6  zero  at  all  hyperbolic  points.  To  solve  this 
system  of  equations  we  have  used  an  Approximate  Factorization  scheme  that  is 
described  in  the  following  section. 

3.4  AFZ  Scheme 

During  each  iteration  of  the  cycle  a  rapid  approximate  factorization  method 
is  used  to  solve  Equation  (8)  with  q'’  fixed.  The  factorization  is  used  in  two- 
dimensional  cylindrical  (constant  r)  surfaces  about  the  rotor  axis.  The  angular 
extent  of  the  grid  on  these  surfaces  depends  on  the  number  of  blades  present.  The 
grid  is  periodic  at  the  two  vertical  edges  and  the  potential  and  q  "  are  also  forced 
to  be  periodic  there. 


32 


The  basic  idea  behind  the  AFZ  scheme  is  to  solve  a  set  of  implicit  equations 
in  each  constant  radius  plane,  for  a  correction  to  4>-  Since  no  implicit  equations 
are  solved  in  radial  direction,  the  three  dimensional  array  of  <i>  values  can  be  stored 
on  disk  and  only  several  6  -  y  planes  of  data  need  be  stored  in  the  computer  at 
any  one  time. 

In  Figure  7,  the  marching  direction  in  the  computational  plane  is  shown. 
The  equations  are  solved  separately  in  each  half  plane  (upper  and  lower)  with 
boundary  conditions  that  ensure  continuity  across  the  interface  at  convergence. 
During  the  iteration  sweep  number  n  +  1,  when  updating  plane  number  k,  values 
of  4>  are  available  corresponding  to  iteration  n  +  1  in  plane  k  -  \  since  it  has  j\ist 
been  updated.  Only  the  updated  (level  n)  values  are  available  in  planes  k  and 
k  -I-  1.  It  was  decided  to  make  use  of  the  available  updated  values  in  calculating 
residual  at  plane  k.  It  is  most  efficient  to  compute  contravariant  velocities  only 
once  each  iteration  between  each  plane  and  use  them  in  computing  Lo{(p)  for  both 
planes  on  either  side.  To  directly  incorporate  updated  values  of  <t>  into  each  Lo(0) 
computation,  it  would  be  necessary  to  abandon  this  approach  and  recompute  these 
velocities,  using  them  only  once  for  each  computation.  This  would  almost  double 
the  number  of  calculations.  An  alternative,  which  was  chosen  instead,  involved 
adding  the  correction  multiplied  by  an  appropriate  constant  to  Lo(<i>),  computed 
using  only  old  values  at  each  plane.  The  final  factorized  scheme  can  be  written  as 

=  aujL{<i>^)  +  aujCE;6<i)'^ 

=  ,i>^ +  (10) 

where  N^y  is  an  operator  in  the  XY  plane,  a  is  an  operator  in  the  supersonic 
zone  and  a  number  in  the  subsonic  zone,  w  is  a  relaxation  factor  and  E~  is  the 
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shifting  operator: 


6<^k  =  6(t>k-i 


N^y  was  chosen  from  a  successful  two  dimensional  ADI  scheme  [30]; 

^xr  =  («x  -  “  KESy)  (ll) 

The  values  of  A  and  B  are  given  by  Equation  (7).  Also,  and  Uy  are  numbers 
at  subsonic  speed,  and  become  operators  for  supersonic  flow. 

To  approximate  first  order  artificial  viscosity  terms 

a.^  =  a  -  S;  ^lU^6l 


Oy  =  a 


For  flow  in  the  +X  direction  (in  this  mapping  there  is  no  flow  in  -X  direction 
in  the  computational  plane); 


where 


a  =  oo  + 


ao  =  mai{/?o,/?i) 

2 

do  =  Pi  AZ  r[Z)  max  (0, 1  - 

a* 

di  =2  Pi  r[z)  C 

4 

ai  =  P4  Pi  AZ  r(Z)  min(l, 
ot2  =  P5  Pi  AZ  r(Z)  m*n(l, 

a’ 

C  =  P2  ao 
r(Z)  =  l-t-rfc/rtip 
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Pi  -  Ps  are  ADI  parameters.  The  relaxation  factor  u  is  given  by 

a;  =  Ps{l  +  l/r(Z)) 

The  form  of  oq  was  chosen  to  become  small  when  the  local  Mach  number  ap¬ 
proached  1,  which  is  required  for  stability.  The  use  of  l/r(Z)  in  the  formula  for  w 
results  in  a  relative  under-relaxation  for  larger  values  of  r,  which  further  increases 
stability  and  improves  convergence. 

The  solution  sequence  consists  of  a  set  of  cycles.  In  each  cycle  there  are 
Jt  sweeps  through  the  field,  in  which  values  of  Pi  are  cycled  and  set  equal  to 
P^,eP°,e^P°  •  ■  where  e  is  a  constant.  All  other  parameters  are  kept 

constant.  Best  results  were  found  for  /f  =  4  and  e  =  |  for  fine  grid  calculations. 
In  each  sweep  the  factorized  equations  are  solved  plane-by-plane  starting  from  the 
root,  for  corrections  which  are  added  tp  <i>.  In  each  plane,  first 

(a^  -  6^A6^)4>  —  au}CE~6<f>^  auiL[(p)^ 

is  solved  row-by-row  for  a  temporary  two-dimensional  variable  <t>.  Then, 

(a^  -  6yB6y)6<i>  =  $ 

is  solved  column-by-column  for  the  corrections  6<p. 

3.5  Vortical  Velocity  Calculation 

The  calculation  of  involves  two  parts;  first  the  vortex  sheet  position  is 
integrated  as  a  set  of  marker  streamlines  to  follow  the  flow  using  interpolated 
values  of  ?  from  the  fixed  grid;  then,  from  these  marker  positions  ?  "  is  computed. 
As  discussed  in  section  2.3.3  in  the  decomposition  of  the  potential  and  vortical 
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components,  the  vortical  component,  g  ",  is  not  unique.  The  basic  idea  is  to 
spread  the  sheet  into  a  thin  vorticity  region  described  by  an  appropriate  velocity 
distribution,  g 

After  solving  for  4>  as  described  in  Section  3.4,  a  velocity 

g  =  V(^  +  g  “ 

is  computed  at  the  grid  points.  The  velocity  is  then  interpolated  to  the  marker 
locations.  This  is  done  by  using  a  bilinear  interpolation  of  the  velocities  at  the 
four  nearest  grid  points  to  each  marker  location  (Figure  8).  The  markers  are 
defined  at  each  tf-like  grid  plane  (approximately  constant  0).  Values  of  r  and  y 
are  computed  for  each  marker  using  a  forward  Euler  integration  scheme  marching 
from  one  d-plane  to  the  next. 

rt,(s =  r„(5)  ^  AS  X  g'^Ty)  (12) 

In  Equation  (12)  i\(5)  is  the  current  marker  position.  ^Fy)  is  the  interpolated 
velocity  at  the  marker  location.  To  generate  the  sequence  of  values  for  ry(5),  it 
is  necessary  to  determine  which  computational  cell  each  vector  is  in  and  compute 
an  interpolated  value  of  ^r„)  using  values  of  q  at  corners  of  that  cell.  A  simple 
test  on  r  is  used  to  determine  the  r-plane,  or  K  index  of  the  cell.  The  other 
indices  required  tests  on  the  cross  products  between  vectors  from  the  corners  of 
the  cell  and  r„  and  vectors  between  adjacent  corners.  For  example,  if  the  current 
ry(5  =  5o)  were  in  cell  (j, k),  a  test  that  it  had  crossed  into  cell  (j  l.A:)  upon 
integration  over  one  step  is  that  the  0  component  of  the  cross  product, 

.C'(^)  =  k  -  k^i)  ^  (>-"v(5)  -  r"*,.  k~\) 

changed  sign.  That  is, 

C(5o),  X  C(So  +  A5),  <  0 
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Finally,  the  step  size  AS  is  determined  using  the  interpolated  value  of  A©  and 
r©  (tangential  velocity)  at  the  marker  locations.  AS  is  given  by 

AS  =  — 

r© 

The  values  of  A©  are  obtained  at  the  grid  points  using  the  values  of  ©  at  the 
(»  +  1)  plane  and  the  i  plane.  Once  r„(S  +  AS)  is  calculated  using  Equation  (12), 
the  solution  is  advanced  to  the  i  +  1,  ©-plane.  The  procedure  is  repeated  until 
the  markers  reach  the  end  of  the  grid,  determined  by  ©max  =  ]^(-©mai  <  0  < 
©max)-  In  hover  the  effects  of  other  blades  could  be  exactly  accounted  for  by 
setting  periodic  boundary  conditions  at  proper  angles  given  by  ©max-  Also  the 
other  spiral  turns  of  the  vortex  sheets  can  be  computed  by  relocating  the  markers 
at  the  upstream  edge,  (-(-©max)  th®  same  radial  and  axial  position  where  they 
leave  the  downstream  edge  (-©jtjox)- 

No  instabilities  were  observed  in  regions  of  large  step  sizes.  Thus,  a  large 
number  of  sheets  can  be  computed  and  this  is  restricted  only  by  the  axial  extent 
of  the  grid  and  computer  memory  related  restrictions.  In  general,  an  azimuthal 
length  of  two  and  a  half  times  the  azimuthal  distance  between  neighboring  blades 
is  used  for  computation.  A  computed  wake  geometry  for  a  two  bladed  rotor  and 
a  four  bladed  rotor  are  shown  in  Figure  9  and  Figure  10,  respectively. 

Once  the  set  of  r„  's  are  computed  the  vortical  velocity,  q  "  can  be  computed, 
q  "  can  take  many  forms,  as  long  as  it  has  the  correct  vorticity.  The  best  type  of  ^ 
for  this  purpose  is  the  one  which  is  concentrated  near  the  line  of  vorticity,  so  that 
corrections  to  the  potential  equation  are  limited  to  a  relatively  small  region.  The 
tangential  velocity  changes  rapidly  and  normal  velocity  is  small  near  the  sheet,  and 
a  q  "  was  initially  sought  with  these  properties.  In  this  type  of  representation  a 
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Figure  10:  Cross  Section  of  Wake  Geometry  for  Four-Bladed  Rotor 
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tangential  velocity  discontinuity  is  being  replaced  by  a  region  of  spread  tangential 
velocity  variation.  While  this  is  a  good  picture  for  the  total  velocity,  q,  it  would 
be  a  bad  choice  for  q".  This  is  due  to  the  fact  that  such  a  tangential  q'^  must 
extend  throughout  the  entire  flow  field.  Even  though  the  J  field  is  nonzero  in  only 
a  small  region,  this  particular  q^  field  (which  defines  J)  is  not  small  nor  limited 
to  a  small  region,  and  requires  a  large  amount  of  memory.  Denoting  the  direction 
along  the  normal  to  the  sheet  by  n  and  along  the  line  by  s, 

J  =  dsQn  -  d„q^  =  d^q^'’  -  (13) 

Since  V<^  does  not  contribute  to  u.  Thus,  if  q^  =  0  and  if  q^  goes  to  zero  away 
from  the  line, 

J  F  0 

Hence,  a  9  "  oriented  along  the  sheet  must  be  extended  throughout  the  field  if  the 
vorticity  is  to  be  concentrated  only  near  the  sheet. 

.\nother.  more  computationally  efficient  approach  is  to  use  a  distribution  of 
spread  velocity  normal  to  the  sheet  surface  tis  shown  in  Figure  11.  If  q"  =  0  and 
7^  0,  however.  Jj  —  d^q^  and  g"  can  be  concentrated  near  the  sheet,  going  to 
zero  away  from  it.  Thus,  instead  of  point  vortices  the  sheet  is  considered  to  be 
made  up  of  a  set  of  discrete  velocity  vectors  oriented  normal  to  the  line  within  the 
cross-plane  and  concentrated  near  it.  Then  q  ''  is  non  zero  only  in  a  small  region 
near  the  sheet.  This  representation  seems  very  unphysical  until  one  considers  that 
it  is  not  the  entire  velocity  which  is  represented,  but  merely  the  vortical  portion, 
g  The  total  velocity,  g  must  still  be  of  the  form  shown  in  Figure  12.  after  the 
potential  part,  V(j>  is  added.  The  vortical  part  g can  have  any  form  as  long  as 
J  =  V  X  g“.  The  normal  form  depicted  in  Figure  11  has  the  advantage  of  being 
nonzero  only  in  a  thin  region,  just  like  cJ. 
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Figure  11:  Vortical  Velocity  Distribution  for  the  third  sheet  of  Figure  10 
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Figure  12:  Velocity  in  a  Cross-Plane 
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The  advantage  of  this  type  of  representation  is  twofold: 

1.  It  is  possible  to  give  a  structure  to  the  vorticity  field  u.  This  field  can  be 
either  solved  for  or  modeled. 

2.  It  is  possible  to  put  vorticity  anywhere  in  the  field  with  no  constraints  imposed 
by  the  grid. 

With  these  features,  a  combined  Eulerian-Lagrangian  code  hats  been  developed 
for  the  hover  problem.  The  vortical  velocity  component,  is  a  normal  velocity 
distribution  describing  a  thin  shed  vorticity  sheet  which  is  allowed  to  convect  freely 
through  the  flow  field,  ^"is  a  continuous  vortical  velocity  field  which  suffers  no 
artificial  numerical  diffusion  and  can  eaisily  be  tailored  to  the  problem.  For  the 
hover  problem  it  has  sufficed  to  model  g ‘'using  functions  which  were  chosen  on 
the  basis  of  numerical  convenience  alone. 

In  order  to  find  the  required  strength  of  g ‘’we  use  Gauss'  theorem  to  obtain 
a  relation  for  the  integral  of  along  a  normal  at  each  point  on  the  sheet; 

I  qldn  =  r 

The  circulation,  F,  is  known  at  the  upstream  edge  of  the  sheet  (blade  trailing 
edge)  from  the  lift  distribution  (which  is  computed  as  part  of  the  entire  calcula¬ 
tion).  Since  it  is  constant  along  mean  streamlines  within  the  sheet  it  can  easily  be 
computed  on  the  entire  sheet.  This  relation  provides  a  scaling  factor  which  gives 
the  magnitude  of  g  "  as  soon  as  the  width  of  the  layer  and  the  functional  form  of 
J  is  determined.  This  width  and  functional  form  can  be  found  by  using  a  viscous 
solution  or  by  simply  choosing  apropriate  values.  For  the  rotor  wake  problem,  the 
latter  approach  suffices.  These  vectors  are  convected  with  the  flow,  and  as  the 


line  evolves  and  the  distance  between  vectors  changes,  the  magnitude  of  the  q 
vectors  (<?")  is  changed  so  that  the  averaged  9"  remains  constant  at  each  position 
on  the  line.  To  avoid  singularities  9”  b  spread  over  a  small  region,  comparable  to 
the  mesh  speicing  on  either  side  of  the  line: 

t 

where  is  the  circulation  at  the  trailing  edge  where  the  particular  marker  (i) 
originated,  the  spreading  function, 

Ar^ 

<r(Ar)  =  max{0, 1 - —) 

and  q'l’  is  a  vector  normal  to  each  panel  and  proportional  to  the  panel  area.  The 
normalization  constant  C  is  such  that,  in  the  limit  when  markers  are  closely  spaced 
and  is  slowly  varying,  the  integral  of  9 through  the  sheet  at  any  marker, 
m,  equals  to  Fm-  The  only  specified  quantity  is  ”0”  which  is  taken  to  be  about 
two  cell  widths.  Stable,  smooth  solutions  have  been  found  using  this  spreading. 

Referring  to  Figure  13,  the  closed  line  integral  oi  q-  di  around  a  portion  of 
the  sheet  is  equal  to  the  area  integral  of  Wn-  Since  V<^  does  not  contribute  to 
the  line  integral  and  q''  is  short  range,  only  those  portions  of  the  line  near  the 
sheet  encircled  by  dashed  curves  contribute  to  the  integral.  'No  matter  how  q  ^ 
is  distributed  along  the  sheet  within  the  closed  curve,  the  area  integral  of  uj„ 
is  determined  by  the  two  small  areas  where  the  path  crosses,  which  are  equal 
to  the  circulation  corresponding  to  the  markers  at  those  points.  The  motion  of 
the  markers  through  the  velocity  field  is  determined  by  the  requirement  that  the 
integral  of  momentum  through  the  sheet  is  conserved.  It  is  important  to  see  that 
only  such  short-range  (nonzero  in  only  a  few  cells  on  either  side  of  the  sheet) 
9  "’s  can  be  effective  for  use  in  this  scheme.  If  we  had  a  long  range  form,  given, 


VftUcity 

(||or*4l  to 
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Figure  13:  Schematic  of  Circulation  Integral 
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for  example,  by  the  Biot-Savart  law  then  the  q"  calculation  would  have  to  be 
made  for  each  point  in  the  entire  grid.  Each  of  these  would  effectively  involve 
a  two  dimensional  integration  over  the  entire  vortex  sheet,  which  would  then  be 
prohibitively  expensive. 

An  interpretation  of  the  q*'  term  for  incompressible  flow  is  as  follows:  con¬ 
sidering  a  small  region  near  the  line,  we  can  ignore  the  d^<i>  tangential  terms  in 
the  Laplacian  since  will  be  much  larger  and  we  have, 

V-(V,^  +  9'')  =  VV  +  dn?"  =0. 

or 

-  dl4>  =  -d^ql 

Hence, 

dn<t>  =  -(In  ^  constant 

'M  =  J 

where  brackets  denote  the  jump  in  0  across  the  line  and  the  integral  is  done  normal 
to  and  through  the  line.  Thus  we  are  effectively  putting  in  a  smeared  discontinuity 
in  <p.  Exactly  the  same  argument  applies  to  the  compressible  potential  flow  case. 

The  vortical  velocity  is  computed  at  each  point  in  the  grid  where.  V<p  is 
computed  using  a  spreading  distance  that  is  proportional  to  the  local  cell  area  for 
each  sheet.  Then,  the  contribution  to  q'’  from  each  sheet  at  a  given  point  is  added 
to  give  the  total  q 


3.5.1  Revised  Formulation  for  Vortical  Velocity  Calculation 

The  calculations  described  above  were  done  with  an  initial  formulation  de¬ 
scribed  in  Section  3.1.  It  can  be  shown  that,  as  the  grid  is  refined,  the  (specified) 
vortex  spreading  with  this  initial  formulation  must  decrease  less  rapidly  than  the 
cell  size.  Otherwise,  numerical  effects  appear  which  can  be  reflected  in  the  motion 
of  the  vortex.  For  example,  considering  a  blade  segment  with  constant  circulation, 
if  the  spreading  is  of  the  order  of  the  cell  size,  h,  and  the  spreading  function  is 
constant  within  h  and  zero  outside. 


\q'>]T/h 

if  a  grid  point  is  within  the  spreading  distance  and  zero  otherwise. 

The  vorticity  in  a  computational  cell  where  only  one  point  is  within  spreading 
distance  of  the  sheet  is  then  given  by 

Where  a  is  the  angle  of  orientation  of  the  sheet  and 

0  =  cos  a  —  sin  a 


From  this, 

]wl  =  dT/2h^ 

or 

ju/j  ~  l/h^ 


and  the  numerical  error  in  velocity  even  in  a  region  of  zero 
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Thus,  a  numerical  vorticity  is  created  locally  (when  it  should  actually  be  zero  for 
constant  F)  that  increases  as  This  also  happens  in  regions  where  three  corners 
of  the  cell  are  within  the  spreading  distance  and  one  outside. 

A  more  accurate  form  for  g  "  based  on  Clebsch-type  of  representation  has 
been  developed  and  used  on  four-blade  rotor  configurations.  In  this  formulation; 

r  =  r=V0, 

where  r‘^(^  is  a  three  dimensional  field  which  smoothly  goes  to  the  appropriate 
r  (circulation  value)  on  the  sheet  as  F  approaches  the  sheet  surface.  The  Clebsch 
potential,  smoothly  goes  from  -i-  1/2  on  one  side  of  the  sheet  to  -  1/2  on  the 
other.  A  convenient  formula  for  0(f)  is 

0(f)  =  1/2  sin  (7rSn/2),  |5„1  <  a/2 
0(r)  =  4-1/2,  5„  >  -+-a/2 
0(r-)  =  -1/2,  Sn  <  -a/2 

where  S„  is  a  (signed)  normal  distance  from  the  point  fto  the  sheet. 


We  use  interpolation-like  formulae  to  compute  r‘^(f)  and  5n(f)  at  any  grid 


point  r; 


I 

^n(^  =  (r,fi)<7(Af<)l//l 

t 

s„(f)  =  !5„(f)| 

A  =  J^a(Af«), 


Afi  =  r  -  ft 

where  the  Fj  is  the  circulation  defined  for  each  marker,  /,  that  defines  the  sheet.  F« 
remains  constant  on  each  marker  from  the  point  on  the  blade  trailing  edge  where 


it  is  shed.  Also,  5”(F,  f^)  is  the  normal  distance  from  f^to  a  plane  through  the 
sheet  at  marker  t  and  the  spreading  function,  <T(Af^)  is 

o{^f^  =  mai(0, 1  -  Afifa^) 
where  a  is  a  specified  spreading  distance. 

Even  through  0  is  non-zero  throughout  the  field  (except  on  the  sheet),  VV’ 
is  zero  beyond  the  spreading  distance  from  the  sheet.  Accordingly,  both  r‘^(^ 
and  need  be  computed  only  on  those  grid  points  that  are  in  a  narrow  band 
about  the  sheet  with  thickness  of  the  order  of  the  spreading  distance.  With  this 
<7  ”  there  is  no  spurious  numerical  vorticity  in  regions  near  the  sheet  where  F  is 
constant  (no  physical  vorticity  on  the  sheet),  and  even  in  regions  where  F  varies 
smoothly  there  is  no  divergent  error  that  becomes  large  as  the  grid  is  refined.  In 
regions  where  F  is  constant  F^(f)  will  be  constant  and  q  "  can  be  written 

9"  =  V(FV(r)). 

Even  though  q  "  is  still  non-zero,  if  the  same  numerical  differencing  scheme  is  used 
for  VV'  as  is  used  for  V<^,  the  effect  of  the  sheet  in  regions  of  zero  vorticity  will  be 
identically  zero.  This  was  not  true  for  the  earlier  formulation  of  g  ",  and  resulted 
in  a  requirement  for  larger  spreading  distances  and  thicker  vortices  (for  a  given 
grid).  The  new  formulation  has  been  successfully  applied  to  two  four-bladed  rotor 
configurations  and  is  seen  to  give  improved  vortex  geometry.  Vorticity  contours 
for  a  four-blade  rotor  obtained  using  this  revised  formulation  is  shown  in  Figure 
14.  Also  this  formulation  allows  the  embedding  of  a  tighter  vortex  sheet  with 
smaller  spreading,  even  in  a  fairly  coarse,  computationally  efficient  mesh. 


3.6  Boundary  Condition 


For  all  of  the  results  presented  a  128  x  40  x  32  grid  computational  was  used. 
The  upper  boundary  is  at  about  a  half  rotor  span  above  the  rotor  disk.  Dirichlet 
conditions  were  imposed  on  the  potential  there,  using  a  form  given  in  Reference 
[14]  for  a  semi- infinite  vortex  cylinder,  which  approximates  the  vortex  system  at 
this  boundary,  given  by 


where  y  is  the  position  of  the  top  of  the  cylinder.  At  the  lower  boundary,  about 
0.8  span  below  the  rotor  disk,  a  corrected  periodic  condition  is  imposed  on  both  q  " 
and  the  potential.  In  this  region  the  markers  are  descending  at  an  approximately 
constant  rate  and  the  vortex  sheet 'is  approximately  periodic  in  height  (y).  If  this 
vortex  system  were  infinite  this  velocity  would  be  exactly  periodic.  A  correction  is 
required  since  the  vortex  system  is  only  semi-infinite,  starting  approximately  from 
the  rotor  plane.  This  correction  involves  adding  a  velocity  similar  to  that  used  at 
the  upper  boundary. 

At  the  blade  surface,  the  normal  velocity  is  zero,  which  implies 

dn<i>=-q:. 

Finally,  Dirichlet  conditions  are  used  at  the  outer  cylindrical  boundary. 
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CHAPTER  4 

RESULTS 


4.1  Comparison  of  Results  for  Two-Blade  Rotors 


Several  applications  have  been  completed  which  illustrate  the  predictive  ca¬ 
pabilities  of  the  method.  Some  of  the  results  reported  for  the  two  bladed  rotors, 
were  published  in  1986  [32,  33,  34].  Results  are  presented  for  the  full  solution  of 
a  lifting  rotor  in  hover,  as  computed  by  the  three  dimensional  compressible  flow 
code,  HELIX  I  which  has  been  developed.  In  Figure  12  preliminary  results  of 
the  full  three  dimensional  compressible  code  eu'e  presented  for  a  single  vortex  sheet 
shed  by  a  rotor  blade.  Here,  the  velocity  in  a  cross-plane,  normal  to  the  blade 
motion,  several  chords  behind  the  trailing  edge  is  presented.  The  cross-stream 
velocity  clearly  shows  a  tight  vortex  which  is  the  size  specified  by  the  model  used. 
This  remains  the  same  as  it  is  convected  downstream.  Other  internal  vortex  struc¬ 
ture  and  density  distributions  consistent  with  the  grid  spacing  could  be  imposed. 
For  this  analysis,  a  rotor  with  an  aspect  ratio  of  19,  a  constant  pitch  of  5*^  and 
Joukowski  profile  was  used.  Figure  15  depicts  the  vorticity  contours  for  the  sepa¬ 
rate  contribution  of  each  individual  sheet  for  relative  values  of  0.45  and  0.7.  The 
contour  values  are  chosen  to  show  the  effect  of  the  individual  sheets  as  well  as  the 
spreading.  Figure  16  presents  the  contours  for  the  total  vorticity  field,  which  is  a 


sum  of  individual  sheet  contributions 


The  next  two  bladed  rotor  chosen  for  analysis  was  a  high  aspect  ratio  {AR  = 
18.2)  rotor  with  a  constant  8®  pitch  and  NACA0012  profile.  Flow  visualization 


Figure  15:  (continued) 
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data  are  available  for  this  rotor  from  Reference  [4j,  and  the  calculated  cases  were 
matched  to  this  reference.  The  computed  wake  geometry  along  with  the  experi¬ 
mental  data  in  the  form  of  tip  vortex  axial  and  radial  coordinate  are  presented  in 
Figure  17.  The  load  distribution  is  presented  in  Figure  18  and  chordwise  pressure 
distribution  at  68%  and  89%  span  are  presented  in  Figure  19. 

A  final  two-blade  rotor  case  has  an  aspect  ratio  of  6,  NACA0012  profile  with 
no  twist  or  taper  a  collective  pitch  of  8°.  The  experimental  data  used  for  compari¬ 
son  was  obtained  from  Reference  [35).  In  this  report  a  complete  set  of  experimental 
data  for  chordwise  pressure  distribution,  tip  vortex  geometry  and  spanwise  load 
distribution  was  provided  for  a  wide  range  of  tip  Mach  numbers  including  tran¬ 
sonic  fiow.  Two  tip  speeds  were  considered  in  this  computation:  one  subsonic 
(Mtip  =  0.436)  and  another  transonic  (0.877). 

For  the  subsonic  case,  the  computed  wake  geometry  is  presented  in  Figure  20 
as  the  tip  vortex  radial  coordinate  and  axial  coordinate  plotted  as  a  function  of 
vortex  age.  The  experimental  vortex  geometry  of  Reference  [34]  is  also  presented  in 
Figure  20.  It  can  be  seen  that  the  computed  vortex  geometry  compares  favorably 
with  experimental  data.  Also  it  can  be  observed  from  Figure  17  and  Figure  20  that 
the  wake  parameter  k  (axial  slope  of  the  tip  vortex  trajectory  before  passage  of 
the  following  blade)  shows  a  strong  dependence  on  aspect  ratio.  This  is  expected 
since  for  rotors  of  the  same  blade  number  and  operating  condition,  both  the  disk 
loading  Ct  and  infiow  velocity  will  decrease  with  increasing  aspect  ratio,  resulting 
in  lower  tip  vortex  settling  rate.  In  Figure  21,  the  circulation  distribution  is 
presented.  There  is  a  relatively  large  increase  in  the  normal  loads  due  to  the  tip 
vortex  shedding  in  the  last  20%  of  the  blade  radius.  This  shows  the  importance 
of  modelling  of  the  tip  vortex  shedding  for  prediction  of  rotor  loads.  Chordwise 


Figure  17:  (continued) 


CIRCULATION  ACROSS  BLADE  <AR-ie.  2> b-2> 
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Figure  21;  Blade  Circulalion  Distribution  for  Low  Aspect  Ratio  Two-Bladed 
Rotor  with  a  Collective  Pitch  of  8'^  and  NACA00I2  Profile  at 
Subsonic  Tip  Speed 
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surface  pressure  distributions  for  the  68%  and  89%  radial  locations  compare  very 
well  with  experimentally  measured  data  (Figure  22). 

The  computed  solutions  for  trsuisonic  tip  speed  are  presented  in  Figures  23 
through  Figure  25.  The  computed  tip  vortex  axial  coordinate  and  radial  con¬ 
traction  are  given  in  Figure  23.  The  bound  circulation  for  this  rotor  is  given  in 
Figure  24.  Chordwise  surface  pressure  distributions  compare  favorably  with  the 
experimental  data  at  all  radial  stations.  Figure  25  illustrates  the  data  correlation 
at  the  68%,  80%,  89%  and  96%  radial  locations.  The  calculated  suction  pressure, 
in  regions  near  the  region  of  tip  vortex  passage  are  accurately  predicted  for  all  but 
the  last  station  which  is  very  close  to  the  tip.  To  resolve  this  region  a  finer  grid  is 
required  there.  From,  Figure  22  and  Figure  25  it  is  easily  seen  that  inboard  pres- 
sure  distributions  are  only  slightly  affected  by  rotor  speed.  However,  the  outboard 
sections  show  corisiderable  pressure  alteration  and  shock  development  as  the  tip 
Mach  number  approaches  sonic  values. 

4.2  Comparison  of  Results  for  Four-Blade  Rotors 

Even  though  good  results  were  obtained  with  two-blade  rotors,  four-blade 
rotors  present  a  crucial  test  for  HELIX  I  since  there  the  vortex  from  each  blade 
passes  much  closer  to  the  preceding  blade.  Also  most  modem  rotors  have  four 
blades.  As  discussed  above,  the  two-blade  rotor  results  presented  were  done  with 
an  initial  formulation  of  the  vortex  embedding  method.  Results  for  four-blade 
rotors  were  obtained  using  the  new  formulation  described  in  Section  3.5.1.  The 
first  four-blade  case  computed  involved  blade  number  7  of  Reference  [35|.  The 
apsect  ratio  of  this  rotor  is  15,  with  a  0A209  profile  and  a  collective  pitch  of  10°. 
The  blade  has  a  linear  twist  such  that  the  pitch  is  increased  by  8.3°  from  root  to 


65 


0.  000  0.  500  1.  000 

X/C 

CHORDW I SE  PRESSURE  DISTRIBUTION  <AR  -8,  b  =2.  r-ZR-B.  68) 


X/C 

CHORDWISE  PRESSURE  DISTRIBUTION  CAR  -8.  b  »2.  r-/R-0,  89) 

Figure  22:  Comparison  of  Chordwise  Pressure  Distribution  for  Low  Aspect 
Ratio  Two-Btaded  Rotor  with  a  Collective  Pitch  of  8*^,  NAC  A0012 
Profile  at  Subsonic  Tip  Speed  with  Experimental  Data  of 
Caradonna  [34] 
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Figure  23:  Comparison  of  Tip  Vortex  Coordinates  for  Low  Aspect  Ratio 
Two-Bladed  Rotor  with  a  Collective  of  S*’  and  NACA0012  Profile 
at  Transonic  Tip  Speed  with  Experimental  Data  of 
Caradonna  [34] 
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Figure  23:  (continued) 
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Figure  25:  Comparison  of  Chordwise  Pressure  Distribution  for  Low  Aspect 
Ratio  Two-Bladed  Rotor  with  a  Collective  Pitch  of  8*^,  NACA0012 
Profile  at  Transonic  Tip  Speed  with  Experimental  Data  of 
Caradonna  [34 1 
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tip.  In  Figure  26  the  computed  vortex  geometry  together  with  experimental  data 
is  presented.  In  Figure  27  the  computed  load  distribution  is  plotted  as  a  function 
of  span  and  compared  with  experiment.  Where  experimental  results  are  available 
agreement  is  similar  to  that  obtained  in  the  two-blade  cases. 

Finally  a  high  aspect  ratio  {AR  —  18.2)  four-blade  rotor  configuration  with 
a  NACA0012  profile,  —8®  linear  twist  and  a  collective  pitch  of  8°  was  analyzed. 
The  experimental  data  used  for  comparison  were  obtained  from  Reference  [4].  As 
the  number  of  blades  increase,  each  blade  is  closer  to  the  tip  vortex  generated  by 
the  blade  ahead  due  to  the  reduced  separation  between  the  blades.  This  causes 
an  increased  aerodynamic  interference  and  a  steep  increase  in  load  distribution 
near  the  tip  region.  This  can  be  clearly  seen  in  Figure  28,  where  the  load  distri¬ 
bution  is  presented  as  a  function  of  normalized  radial  distance.  The  computed  tip 
vortex  axial  and  radial  coordinates  are  presented  in  Figure  29.  Agreement  with 
experiment  is  seen  to  be  good  in  this  case. 

4.3  General  Features  of  the  Solution  and  Wake  Geometry 

Several  general  features  of  the  tip  vortex  geometry  are  evident  in  Figures  20 
through  29.  When  an  element  of  the  tip  vortex  is  shed  from  a  blade,  its  axial  rate 
of  descent  is  low  until  it  passes  beneath  the  following  blade.  At  that  point,  the  tip 
vortex  element  lies  radially  inboard  of  the  tip  vortex  of  the  following  blade  and  thus 
experiences  a  large  downward  induced  velocity.  The  axial  tranport  velocities  before 
and  after  the  passage  of  the  following  blade  are  fairiy  constant  in  the  near  wake, 
«is  can  be  seen  from  the  substantially  linear  variations  of  the  axial  displacement, 

with  wake  azimuth  angle  in  these  regions.  The  radial  displacement,  r^,  of  the 
tip  vortex  decays  in  approximately  exponentially  with  increasing  azimuth  angle. 


I 

1 
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I-  igure  26:  Comparison  of  Tip  Vortex  Coordinates  for  a  Four-Bladed  High 
Aspect  RaUo  Rotor  with  8.3*^  Linear  Twist  and  10*^  Collective 
rjtch  and  OA209  Profile  with  Experimental  Data  of  Maresca  [35j 
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Figure  27:  Comparison  of  Blade  Circulation  Distribution  for  a  Four-Bladed 
High  Aspect  Ratio  Rotor  with  -8.3"  Linear  Twist  and  10"  Collec¬ 
tive  Pitch  and  OA209  Profile  with  Experimental  Data  of  Maresca 
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Figure  28:  Blade  Circulation  Distribution  for  a  Four-Bladed  High  Aspect 
Ratio  -8'^  Rotor  with  Linear  Twist  and  8'^  Collective  Pitch  with 
N  AC  AGO  12  Profile 


Figure  29:  Comparison  of  Tip  Vortex  Coordinates  for  a  Four-Bladed  High 

Aspect  Ratio  Rotor  with  -8*^  Linear  Twist  and  8®  Collective  Pitch  with 
NACA0012  Profile  with  Experimental  Data  of  Landgrebe  [4j 
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The  calculations  presented  were  done  using  HELIX  I  on  a  128  x  42  x  32 
grid  with  80  points  on  the  airfoil  and  18  span  stations  on  the  blade.  Typical 
computation  times  on  a  VAX-11/785  system  were  about  5-6  minutes  per  iteration. 
The  code  has  also  been  run  on  the  NASA  Ames  Advainced  Computational  Facility’s 
CRAY  X-MP  system.  The  CPU  time  for  one  iteration  using  this  system  is  of  the 
order  of  10  secs.  No  special  effort  has  yet  been  made  to  vectorize  the  program  for 
use  on  this  machine.  Substantial  improvements  in  the  required  computation  time 
could  be  realized  by  restructuring  the  code  to  benefit  fully  from  the  CRAY  vector 
processing  capabilities.  It  is  estimated  that  by  doing  this  the  execution  time  could 
be  reduced  by  a  factor  of  3  to  4 

Generally,  a  fully  converged  solution  and  a  stable  wake  geometry  are  obtained 
in  about  3  to  4  wake  and  vortical  velocity  updates,  in  betvyeen  which  an  approx¬ 
imately  converged  potential  solution  is  obtained.  It  is  estimated  that  generation 
of  one  complete  converged  solution  for  loading  and  wake  geometry  takes  about  I  ^ 
hours  of  CPU  time  on  the  CRAY  X-MP. 
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CHAPTER  5 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  main  objectives  of  the  research  were  to  develop  a  comprehensive  method 
for  the  prediction  of  hovering  rotor  performance  in  compressible  flow.  These  have 
been  met.  The  technique  includes  the  all-important  wake  effect  and  allows  for  the 
analysis  of  advanced  rotor  configurations.  As  such,  the  constraints  imposed  by 
earlier  techniques  have  been  removed.  Example  calculations  were  presented  which 
demonstrate  the  capabilities  of  the  method  including  a  stable  calculation  for  the 
free  wake  flow.  This  has  been  demonstrated  for  both  a  traditional  high  aspect 
ratio  blade  and  for  an  unconventional  blade. 

Results  for  the  surface  pressure  distribution,  the  blade  loading  and  the  near 
wake  geometry  are  all  in  good  agreement  with  data  but  the  far  wake  geometry 
should  be  the  subject  of  a  future  study  where  improved  calculations  for  ^  can  be 
incorporated. 

The  basic  algorithm  has  been  demonstrated  as  complete  in  every  respect 
except  for  the  recommended  refinement  of  selected  components.  In  addition,  the 
following  are  recommended: 

1.  Studies  should  be  made  for  rotors  with  non-rectangular  tip  geometries  to 
establish  the  significance  of  this  geometry  on  the  entire  flow  field  calculation. 

2.  The  extension  to  rotors  in  forward  flight  should  be  made.  The  major  mod¬ 
ification  here  is  to  the  potential  flow  calculation  which  must  then  be  based 
upon  the  unsteady  flow  equations. 
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